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MEAN GRADIENT HEAT TRANSFER 
IN ISOTROPIC TURBULENCE 

James C. Hill 


ABSTRACT 

The problem of heat transfer of a constant -property fluid in 
isotropic turbulence in the presence of a uniform mean temperature 
gradient is examined. Three closures of the statistical moment 
equations are discussed: the first-iteration approximation, the third- 
eumulant discard, and the direct-interaction approximation. Evalu- 
ation of the thermal eddy diffusivity for the direct-interaction ap- 
proximation requires solution of a nonlinear integro-differential 
equation for the averaged Green's function for the temperature 
field. The numerical approximation of was attempted in spatial 
coordinates to gain experience for shear-flow problems and so far 
is unsuccessful because of numerical difficulties associated with the 
singular nature of <G]> at zero time. The calculation can best be 
done in wavenumber coordinates in which the Fourier transform of 
<(G)> is initially unity. 
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MEAN GRADIENT HEAT TRANSFER 
IN ISOTROPIC TURBULENCE 


X. INTRODUCTION 

The problem treated here* is the calculation of the steady-state mean 
heat transfer rate for a specified uniform mean temperature gradient super- 
imposed on a steady isotropic turbulent velocity field with known low-order 
statistics (e.g., space-time velocity covariance). The fluid is incompressible 
and has constant physical and thermal properties. The temperature field is 
passive, and viscous dissipation is negligible. 

We denote statistical averages (ensemble averages) by angular brackets, 
so that the temperature and velocity fields for any realization are written 
T(x,t) = <T(x,t)> + # ( x» t) and U(x,t) = <U(x,t))> + u(x,t) where $ and u are 
the fluctuating or random parts of the fields. The heat equation for <(T> reduces 
to <u0> = constant vector, whose value is some function of V<(T)> and of the 
statistics of the velocity field. For the present work we wish to determine the 
eddy diffusivity tensor k 1 , defined by = -K t ’ V <t)> , by making statistical 
or dynamical assumptions about the interaction between the velocity and tem- 
perature fields. In this report several such assumptions are discussed, the most 
interesting to us being Kraichnan f s direct-interaction approximation/ 1, 2) 3 > 
which has not yet been studied for turbulent transport of a passive scalar by a 
mean gradient. Also, the problem is to be done in physical coordinates, with 
the idea of generalizing numerical techniques to shear-flow problems. t 

Although this problem is rather idealized, data are available for a grid- 
generated turbulent flow heated to a self- maintained uniform temperature 
gradient.^ 4 ) 


"Suggested by Prof. C. A. Sleicher. 

'(The theories treated here, including the direct-interaction approximation, are set up in terms of 
Green’s functions which have singular behavior at zero time in the spatial representation. In the 
wave-number or Fourier representation the singularity does not occur but, unfortunately, for non- 
homogeneous turbulence the Fourier representation requires a larger number of convolution 
integrals. 
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2. THEORY 


For the general heat transfer pi’oblem (transient, non-homogeneous) let 
£ be the operator 


c= J_ + <u> • v-«v s 

dt 


where k is the molecular thermal diffusivity. Then the instantaneous temper- 
ature field satisfies £T(x,t) = -u (x,t) *VT(x,t) subject to suitable initial and 
boundary conditions. In terms of the Green’s function^ G(x,t |x', t') for re- 
sponse of the temperature field at (x,t) to delta-function initial conditions, 

G - 8 (x-x') at t = t' , and zero boundary conditions, G = 0 for x or x' t: B , the 
solution for arbitrary initial conditions T 0 (x):x eR and Dirichlet boundary 
conditions T B (x,t):x elB is 


T(x, t) 



x , ' t 0 )T 0 (x')d 3 x'- 


k[ f T B (x', t') -5- G(x, t I x',t')dt'd 2 x' 

hi t 0 Bn 


where n is the outward normal unit vector on 13 . A similar treatment can be 
made for Neumann or Churchill boundary conditions. We shall use the above 
equations in a more symbolic notation with arguments 0=^ (x,t), 1 =£ (x' ,t' ), 

2=^ (x",t"), etc. If the arguments are omitted, the corresponding primes are 
retained (i.e,, u(x",t") - u(2) - u'') unless confusion can arise or if only a 
symbolic representation is intended. The initial and boundary conditions are 
written simply as <I>, so the solution T(x,t), which was written out above in de- 
tail, is simply T = /G(01) <I» (1) d 1. Both T and G are random functions as a 
result of the random coefficient u(x, t). The random fluctuations in <1> and G will 
be denoted by <f> and g . 


The temperature -field equation can be used to formulate expressions for 
statistical moments of T or G . Unfortunately, each equation contains a higher- 
order moment because of the random coefficient u , and an infinite set of un- 
closed equations or an infinite iteration expansion results. The task of devising 
an approximation to form a determinate mathematical problem by suitably 
truncating or consolidating the set of equations or the expansion is called the 
closure problem. 

The methods of closure described below will not be restricted to steady 
isotropic turbulence with uniform V<^T> until Section (2c) where k* is discussed 
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and where the equations are reduced to the isotropic case. In order to define k* 
suitably, it will be necessary to assume that the smallest dimensions of the 
region H is large compared to the scale of turbulence \ ,* so that, over the re- 
gion where we wish to calculate k 1 , V <[T)> is effectively uniform. 

I he Green’s function for the case of zero turbulence is written with a 
superscript zero (since it is a zeroth approximation of G ) and is not a stochastic 
quantity. G° is used below in iteration expansions. For the isotropic case G° 
is the Green’s function for heat conduction in an infinite region 


G°(x,t | x\ t') s [47 ik (t - t')3‘ 3/2 exp {- | x ~x' | 2 /4/<(t -t') > , t> t'. 


a. The maximal randomness condition (MRC) 

The analog of Kraichnan’s MRC* for the convection problem is that the 
statistical dependence between u and T (or G ) is induced solely by the random 
coefficient u and not at all by the boundary or initial conditions. It is sufficient 
to require that u and 0 be uncorrelated to all orders of u . If u and 0 were 
correlated/ then a statistical dependence between 0 and G would be developed 
through u , The MRC prevents the transmission of any statistical information 
of 0 to G . Hence the Green’s function problem can be solved independently of 
the temperature -field problem. We have, then, « 

<T> = j<G(01)> <1>(1)> dl 
< UT > =/< uG(01)><<I>(l)> dl 

<T 2 > = JJ<G(01)G(02)> <t»(2)> dld2 ♦ 


*ln isotropic turbulence define the scale X to be X =/ C ° f(r,o)dr where f(r,t) is the longitudinal 
space-time double velocity correlation coefficient. 0 

^The MRC was originally defined for the Fourier representation of the velocity field in homogeneous 
turbulence.^ The initial conditions must be ''maximally random”; i.e,, the velocity field initially 
is independent, multivariate Gaussian (which follows from central-limit reasoning). A more formal 
statement of the MRC and its use to develop a hierarchy of equations for single-time moments is 

given by Orszag.^ 6 '?) 

tThis would be the case if, for example, part of the boundary IB were an imaginary surface within 
the fluid where there were both velocity and temperature fluctuations. This condition could arise 
for heat transfer to art open system in which we introduce an imaginary boundary to restrict the 
region for numerical computation. The MRC can be satisfied by the trivial case 0 = 0 (i.e., T 0 = 
<T 0 > and T b =<T b ». 
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Tlius the closure problem for the temperature field is equivalent to the closure 
problem for the Green's function if the MRC is specified, We will refer to the 
problem of calculating these temperature -field moments as the T -problem when 
the closure approximations are applied directly t ) the temperature-field and as 
the G-problem when the approximations are applied to the Green’s function. 

In the present problem the MKC need not be explicitly stated since for large 
times and far from the boundaries <u«A(l)>, etc,, 0. We also require that the 
velocity field be statistically homogeneous in space and time and that kd 3 x. 

H 


b. Closures for T and G 


In this section we describe the closures to be compared. The approxima- 
tions used here are the first-iteration approximation (FIA), the third -cumulant 
discard (3CD), and the direct-interaction approximation (DIA), The idea of these 
closures is to separate the statistics of the velocity and temperature fields, so 
that given <uu'> , for example, we can calculate <T> , <nT>, and <T 2 >. These 
closures are of low enough order that (uu') is sufficient information.* Note 
that<^uu'> is an incomplete and non-unique statistical specification of the 
velocity field. Since FIA and 3CD use no variational principle or other maxi- 
mizing device in relating higher to lower moments (only a statistical hypothesis 
on the relation between moments) the higher order statistics of the velocity field 
do not affect the temperature field. Although the DIA is based on a dynamical 
rather than a purely statistical hypothesis and includes no ad hoc assumption 
about the relation between moments, it too results in a temperature field that is 
independent of higher order velocity moments. 

The FIA and 3 CD specify the form of the triple moments containing 0 or g 
as random variables. If the moment is linear in 6 or g, as in the approximation 
for<(u0)> , then the closures for the T-problem are identical to those for the G~ 
problem. The reason is that <l> appears linearly in the G-problem, so 0 averages 
out. However, if the moment is bilinear in 6 or g, as in the approximation for 
6 2 )>, the closures for the T-problem are not the same as those for the G-prob- 
lem. Since appears bilinearly in the G-problem, moments of the form 
appear that do not arise the T-problem. The reason that the T and G problems 
are not identical is because T has stochastic properties (l> initially and on the 
boundary, whereas G does not. Therefore, specification of the form of triple 
moments involving g is statistically less restrictive than specification of moments 


*This selection of closure approximations is chosen to parallel that of Kraicbnan^ in a discussion 
of the asymptotic short and long-time behaviors of random convection in homogeneous turbulence 
for large P^cl^t number. 
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involving 0 . In other words, 0-moments arc intimately tied up with statistical 
properties oM> near the boundaries ofh® (t 0 < t <go), whereas g-moments arc 
not. The exception to this ease, as noted above, is when 6 or g appears linearly, 
and tiie two problems are thus equivalent. In the present problem, for large 
times and far from the boundaries, terms involving go to zero, so that the 
results from the T and 0 problems are identical. 

The invariance and realizability properties of these closure approximations 
have been remarked on by Kraichnan in a different context.^ 8 ) They will not be 
repeated or extended to the more general case described here. 

(1) The FI A. This approximation is equivalent to expanding the temperature- 
field in powers of PecWt number and retaining only the first order term, which 
corresponds to a simple perturbation expansion. The temperature fluctuations 
are expanded about the zeroth order fluctuations by iterating 


0 = 0 ° + J G°(01)V' • [<u'0'> -n'O' - u'<T'> ]dl 


and x'etaining only first order terms. A similar procedure is used for g . An 
alternative approach is to use the following statistical postulate for calculating 
<u0>from the above equation: <uu'0'> = (uu') (T'°- <^T')> ), where T° is 
calculated from G° and Equivalently, <Wg(12)> = <uu'>(G°(12) - <G<12)>). 

The result for <(uS)> In either case is 


<u0>=- Jg 0 (01)< uu'} * V'T'° dl. 


For calculating the needed FIA postulate for the T-problem is <V0'd}= 
f <(u' 0} (T'° - <T'>) where <^u'0} is also evaluated with the FIA. The result is 


<^ 2 >=JrG 0 (Ol)G 0 (O2)<y'u'> : V'T' 0 V"T" 0 dl d2. 


The corresponding G-problem postulate is <V'g (01)g(32)> = <ii”'g ( 01 )>(G° (32) - 
<(G(32)^>) where <^u'"g'01))> i© also evaluated with the FIA. The following term 
must be added to the above result for : 
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JJ|<G(Ol)><G(02)>+JjG 0 (O3)G 0 C04)9'"V IV :<u I V'>G 0 (32)G 0 (41)d3d4|<<Kl)/(2)>c!l t l2 
where <(G> is given by the FIA by 


<G(01)),sG 0 (01) +JjG°(02)<yV')>:7V(23)V"'C? ) (31)d2d3. 


(2) The 3CD. The statistical postulate for this approximation is that all 
triple moments involving $ or g are zero. That is, moments of the form 
<uue?> - <uug^> = <^u 60y - <ugg)>~ 0 for arbitrary arguments of 6, u , and g , 
The results for the T -problem are 


<u 0} = - f G° (01) <uu'> *V'<T'>dl 

<f P} =JJg°(01)G°(02)< u "u'> :V'<T'>V"<T”>dld2 
The G-problem results are identical except ioi the additional term in 


|Jj<G(0l)^G(02)> +JjG°(O3)G°(O4)V , "7 IV :<;u I V , )><(G(32))><^(41)>d3d4|<0(l)^<2))>dld2 


where <^G^ is given by 


<G(01)>=:G 0 (01) +JJg°(02)<u'V'>: 


V w G°( f 23)V" , <G(31)>d2d3. 


(3) The DIA. This approximation has been discussed in detail by Kraich- 
nan / 1 ’ 2 ’ 3) No explicit statistical postulate is made about the relations between 
moments, contrary to the FIA and 3CD, because the approximation concerns the 
dynamics of the random convection rather than the form of a particular statistical 
moment. From Equation 7.3 of Reference (3) the DIA result for<^u0J> is 

= -J <^G(01)]Xiiu'^>* V'^T'^dl 


6 



where <G> is given by 


* 


<G(01))> G°(01) i [Jg°(02)<u^i"^: V''<G(23))>V"'<G(31)>d2(l. 


The calculation of is much more difficult than for the FIA and 3CD, and 
we have not succeeded in developing a closed expression for it. Given <^G)> , 
(uu') t and<(T,> » however, one can calculate <00'^ from the following general- 
ization of Equation 11.33 of Reference (2); 


£<««'):V(t) •J<G(12)> <uu n ) • V''<T">d2+V- /< UU "> 
<G(12)>V ,, <^^ , ^> + <(3(O2)>V"<0 , 0 / '>Jd2. 


Then <0 i 2 )> can be recovered by taking (x' , t')=^ (x,t). Because of the difficulty 
of this problem, the comparison of <^0 2 > for the various approximations will not 
be discussed. This means, unfortunately, that we will be unable to discuss the 
correlation coefficient 


0 = <u0> 


i <u • u><6> 2 > 


- 1/2 


c. Specialization to the present problem. 

Now we introduce the conditions of the present problem. Assume a stationary, 
homogeneous turbulence with <^U^> - 0. Two parallel planes within this field are 
separated by L >> K and are each maintained at a different uniform temperature 
such that the MRC is satisfied (e.g., 0 B & 0). Then for steady state (t - 00) and 
far from the boundaries we expect on physical grounds that a uniform V<T>will 
develop, V <(T^> = A. This, we hope, leads to the existence of an eddy diffusivity 
K t independent of A. The result for <u0)> given by the FIA includes the factor 
VT°= A 0 . The relation between A and A 0 consistent with the FIA is 

A = a - A 0 

where a = 1 + ff VG° (01)<u"u')>- V' G° (12) dld2. For isotropic turbulence a is 
defined by one scalar a , so A = a A 0 . Let us define 


7 




K 0 = j* G°( 01) <(uu'^> dl. 

Then the FIA result is k£ ia - v^/a, and the 3CD result is simply k* cd ~ K . 
Note that<(G^ is not needed for the 3 CD or FIA results. The DIA result is 


k dia s J<p(°l)^<uu'>dl 


where the equation for <^G)> may be rewritten as the solution of 


£ <G(01))>= J<Vu> : V<G(02)>V''<(G(21)> d2sH(01) 


which is the form we shall use for numerical calculation. 

For the case of steady, homogeneous, incompressible, isotropic turbulence, 
<G> reduces to the special form <G (x, t | x', t ')> = <(G( |x-x' | , t-t'))> and 
similarly for G° .* The velocity covariance reduces to the form 


<(u(x, t) u(x', t')y = Q(x -x' , t - t') 


where 


Q(r, t) = |u| 2 rr 


r 

2 


(r,t) 


+ u 


f(r,t) + 


lIL 

2 3r 


(r, t) . 


*This is only true when or G° is multiplied by ']> with the same arguments, so that far 
from the boundaries the spherical approximation is good over distances the order of (Keep in 
mind that <^G> and G° are zero on 13 .) The spherical approximation cannot be used in determining 
a for the FIA, since in the integral expression for a there is one G° not weighted by ^uu'^ 
with corresponding arguments. In that case we must use the nonhomogeneous G° for parallel 
planes separated by L » Let e 1 be the direction normal to the planes with coordinate x, r = 
|(x - x ' ) x e 1 | , and t > t ' . Then G° takes the form 

00 

„ o / i . s 1 V - ’ . /rmxN . fnrrx 

G°(x,x ,r, t-t ) sxn(_) Sln 

n = 1 


)1 


d^J 0 (fr)exp 


I"* 


rmi 

hi 
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In this equation |u| 2 = 1/3 <u • u)> , r = |r], r = r/r, and f is the longitudinal 
space-time velocity correlation coefficient. The eddy diffusivity goes to the 
diagonal form k t = 1 k* . For large times ( t -* a> ) k* may be written 


K* =47 T U 




r 2 dr G( r, t) 




where G(r,t) = G°(r,t)/a for the FIA, G(r,t) = G°(r,t) for the 3CD, and G(r,t) = 
<(G(r,t))> for the DIA. Here a will depend on position and the separation distance 
L, indicating that the FIA is not a local enough approximation to given uniform 
k* . <(G( r, t))> for the DIA is the solution of 


f_i _ 1 ± 



<G(r, t)> = H(r, t) 


(1) 


where <(G ( r, 0)^> = S(r)/277 r 2 , <(G(r,t)> ►O, and H is given by 

r—oo 

J wt r 00 fr 2 > ^ 

dt' 2wr'2dr' d (U <G(r',t')>j f(r',t') + ZJ. f(r',t') 
o i v.L 


2L<G(r",t")> + -l. 


where y 2 - (1 - fj?) r 2 /r" 2 ,r" 2 =r 2 + r' 2 - 2/ir r' , t" = t - t', and use has been 
made of spatial homogeneity. The geometry of the space convolution is shown 
in the following figure. 
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d. Input information. 


The calculation of k* requires specification of k , |u | , and f (r, t). Data are 
available for f < r, t ) for the case of nearly isotropic grid -generated turbulence 
and are summarized by Favre.< 9 > The behavior of f is too complex to permit 
efficient calculation of H , so for the present work simplified approximations will 
be used. One simple model that was briefly considered is the M on-off" function 


f t (r t) s 


1, r < X. and t < r 


0 r > \ or t > r 


where K and r are the integral length and time scales. Unfortunately, f t is not 
realizable because the corresponding energy spectrum function 


E(k) 



( 


sin k k 
k k 



- cos kk 


is not positive for all k. Another model is the exponential function f 2 (r,t) » 
exp {-r/A - t/r} , which has an undefined Taylor microscale. More complex and 
realistic models can be introduced if the numerical procedures prove feasible. 


3. NUMERICAL CALCULATION OF<G> 

One can now calculate k t for the FIA and 3CD since G° and f (r,t) are 
given (and a can be evaluated from these). However, the DIA result for k t in- 
volves <[gX which was given above by the implicit equation £ <G)> = H where H 
is a nonlinear functional of <G>. This equation, Equation (1), must be solved by 
numerical approximation. 

We follow the usual method of approximating ^G^ on a space-time mesh as 
the solution of a set of difference equations which are solved by proceeding from 
one time layer to the next.* The set is consistent with Equation (1) as the space- 
time mesh is allowed to become arbitrarily dense. Unfortunately, little more 


*The numerical solution of parabolic partial differential equations is discussed in (10), (11, Ch. 
7), (12, CH. 8), (13, Ch. 2), and (14), in which many references to the literature may be found. 
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can be said about the validity of the approximation. The existence of a unique 
solution to Equation (1) has not been demonstrated, nor have the stability and 
convergence of numerical schemes been investigated. These schemes, similar 
to those for solution of the spherical heat equation with a source, would be known 
to be stable and convergent if H and initial <(G^> were bounded (or else if certain 
derivatives of <fo^> have singularities of sufficiently low order as t - 0), which 
is not the case. Also, it is difficult to reproduce numerically the strong initial 
singularity in<tf>, and hence an infinitesimal-discretization-error analysis may 
not be appropriate for approximations to Equation (1). Nevertheless, we will 
consider a posteriori evidence of convergence and stability plus the consistency 
requirement to be sufficient justification for the suitability of a numerical method. 


It is convenient to dimensionally normalize both the difference and differ- 
ential equations. The length scale is X, so then normalized <(G^> is X 3 <(G^> and 
normalized r is r/X . If t is normalized by viscous decay time, i.e., as tv-'k 2 , 
then in Equation (1) we take k -♦ (Prandtl No.) and |u| -♦ Reynolds No. ( | u | X/v) . 
If t is normalized by conduction time, i.e., as t k/X 2 , then k 1 and | u | -Pe'cle't 
No. (|u|X//<), Or finally, if t is normalized by eddy circulation time, i.e., as 
t | u |/X, then k - (Pdcldt No.)- 1 and |u | - 1. We leave the choice of time scale 
arbitrary. 


a. Difference equations for <^G 

The solution <(G (r n ,t m )^> of Equation (1) is to be approximated by G m (or 
simply G ) for all points (n,m) on the space-time mesh. We take uniform space 
increments h so that r n = nh: n = 0, 1, 2,. . .N. The range of r is finite with 
r N » X such that integrals with respect to r can be truncated at r N with error 
0(h 2 ). The time increments are arbitrary, l - t - t m = 1, 2, . . . M 

' ' ^ m m m- 1 

where = 0. The assignment of suggested by Douglas O0) is ^ m + 1 r i m ~ 
constant b £ 1 (m > 0) which takes into account the expected increase in smooth- 
ness of <(G^> with respect to r as time progresses. The time increments may 
also be written 4 m + l = l x + (1-b) or ^ m+1 = b m for m > 0. 


Here we list some definitions of difference operations on space functions 4> 
and time functions \ p : 


Laplacian (spherical) central difference: 

(1/h 2 ) 


l4U n + 1 


, n > 0 


S 2 0 

n ”n ’ 


(6/h 2 ) - <£ Q ] 


, n = 0 
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Gradient central difference: 


(1/ 2h) W> n + 1 - 0 nml ] , n > 0 


S <b 

n mi 


< 


l(l/2h) [4^ -<t> 2 - 30 o ] , n = 0 


Spatial interpolation: 

& = € <t> n + (1 - £) <£ n + ! . (n = mod (x) , £ = x - n) 


Two-layer average: 


s6 ^i +(i-fl)^», (0<e<i, t 9 = t ra + ^ m+1 ) 


Three -layer average: 


< 3> ^m = | (V'O '/ /m + 1 + |^ m + | (^m + l^m) 'A'"*' *• (^m = f + 1 ] 

Two-layer time-difference: 


A< ( 2 > 1 //” + e = [^"> + ‘ - ^”]/X m + 1 


Three-layer time -difference: 

A< t 3 > 0 m = [(l/b) t/r m + l + (b - 1/b) i// m - b\p n - l ]/2\. 


The errors associated with applying these approximations (see Appendix) are 
all second order with respect to mesh spacing if 0 = 1/2. If 6 4 1/2 then the 
two-layer time difference has a first order truncation error. 

The quadrature formula for H is based on the trapezoidal rule. Similar 
formulae may be obtained for Simpson’s rule, but the error terms will involve 
higher order derivatives. The integrand of H , except for the factor 2rrr' 2 , is 
discretized by approximation of the derivative terms to obtain P(r n , r n ,, /x. , t m , 
t , | G ) defined by 

m 1 ' 


p ( . . . 


* 


G) = G”, 



+ F, 


IS* 


g;-*']} , 


(x 2 ss n 2 + n' 2 - 2n n' /x.) 


where 






f ( r n' ’ 





(1 -7 2 ) 



f <V ■ 



r rhe integration over r' is performed by a modified trapezoidal rule obtained by 
integrating the first order Lagrangian interpolation polynomial times the factor 
2 7T r ' 2 over each interval such that 

* 

2wh3 ( n 2 + lV"' 


Since the time intervals increase monotonic ally, it is desirable to perform the 
time convolution by summing simultaneously from times 0 and t m to l/2 t m H 
t where m 2 - log b [1/2(1 + b m )].* The indicated integration is 

t: m m 2~ 1 

( «A(t f )dt # lp 0 + 2^ <// m2 +'t 1 iftj + (0m' + '/'m-m')- 

J ° m'.l 

The triple sum will be denoted by 2 ni J’ . . Then the right hand side of Equation 

n j m 

(1) is approximated by 


*Since m 2 is not generally an integer, a small correction for the interval (m 2 ~b m 2 +D must b® 
made in practice. 


f r N N- 1 

J 27rr' 2 0(r')dr'% +7rh 3 ^N 2 - —N + ^ 


The /x -integration is over intervals ( A/x) = 2/J , 


c i J-i 

0(/x) d/X% (0 O + 0j)/J + Y r '/'j 

•J-l J 7T1 
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H S = P On' r n ( ' h ' V *»< I G ) 

v? jm' 

which has a truncation error maintained atO(h 2 ) +0(^ 2 ) with the restrictions 
outlined in the Appendix. Unfortunately, even if we succeed in restricting the 
truncation error to this order, the overall discretization error may not be small 
since it obeys a nonlinear integral equation involving <^G)> . 

We consider two types of difference equations for G, (I) The Crank-Nicolson 
form , which is set up by approximating Equation (1) on the fictitious time layer 
m+0, 0^0-1, gives an implicit equation for values of G on time layer m + 1: 

M2) G m + 0 _ K g2 <2)Qin + 0 _ ^ + 6 (2) 

t n n n ' ' 

The truncation error associated with (2) isO(h 2 ) +0(^ m + 1 ) unless 6 - 1/2, for 
which it is 0(h 2 ) )• The latter case is preferred unless stability con- 

siderations require 6 > 1/2.* The mesh points of G appearing in (2) are shown 
schematically in the following figure: those associated with £ are marked as an 
asterisk, and those associated with H are circled. The time increments are 
shown to be uniform for simplicity. 



(II) The multilayer form is obtained by applying tbo discretized heat operator 
to three adjacent time layers, 

A< t 3 )G^ -kS 2 < 3 >Gj; =H^. (3) 


*|f H =0 (no turbulence) Equation (2) is stable and convergent for e > 1 / 2 . 
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The truncation error in this case is 0(h 2 ) +0([£*3 2 )> ’ ('Vm+i)'** The mesh 

points for Equation (3) are laid out below. 
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b. Solution of the difference equations. 

Let g be the vector of values of G for all n on the (m+1) time layer. We 
wish to solve (2) and (3) for g in terms of G for all previous time layers. 
Equations (2) and (3) may be written as 


Bg = d + z {g} (4) 

where the left side represents the discretized part of C <(G)>on time layer (m + 1), 
the vector d comes from the contributions to £ <(G)> and H of all time layers < m , 
and the vector functional z is from the portion of H depending on the (m+1) time 
layer. The matinx B is tri-diagonal and may be inverted directly by the Thomas 
algorithm.* 

If (4) is of the Crank-Nicolson form, the term z is treated as a perturbation. 
A sequence of iterations is set up with the k th iterate g< k > given by 

g( k ) = cjB"' 1 ^d + z{g^ k “ 1 ^}'j + (1 - oS) g( k “0 


where eo is a constant called the relaxation factor. For co > 1 this is the method 
of successive overrelaxation. The zeroth interate is calculated from (4) with the 


*The Thomas algorithm p ‘ is a normalized form of Gaussian elimination (see Appendix). 
It is stable with respect to round-off errors if B is tri-diagonal and diagonally dominant, as it is 

for Equations (2) and (3). 
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argument g 55 G m+l of z replaced by G m . The iteration is corwmued until the 
parameter {max n ||g< k > |-| g< k - l >j|/max n |g< k H} is arbitrarily small, usually 
10 “ 4 . 


If (4) is of the multilayer form, then * n 0 and (4) is solved directly. Equa- 
tion (3), however, applies to three adjacent time layers, and so Equation (2) 
must be used to calculate G 1 , 

The outer boundary condition for G Is that GJj = 0* for all m. The values 
of G for n ~ 0, the inner boundary, are given by (2) or (3). An alternative and 
simpler scheme for G$ would be to take S 0 Cg » 0. 

The initial condition of is supposed to be the spherical Dirac delta- 
function, The numerical approximation is to take G° to be a spherical Gaussian 
with initial variance cr 2 « \ 2 . This is only possible if cr Q is the same order as 
h , and so we expect a large interaction between o* 0 and the truncation error. A 
good approximation to <(g)> can be recovered only by taking cr Q -* 0 where h -* 0 
faster, The practicality of this procedure is still doubtful, 

c. Numerical difficulties. 

The numerical approximation of <Jg)> is subject to the usual difficulties 
experienced with parabolic systems regarding the effect of mesh size on con- 
vergence and stability. In addition, the singularity of the initial condition and 
the nonlinear convolution make error and stability analyses difficult (although 
energy methods may be useful). Small truncation error does not guarantee a 
small total discretization error, even for a stable numerical scheme. Below 
we list details of some expected numerical difficulties associated with the singu- 
lar nature of gS . 

(1) The truncation errors for (2) and (3) have terms 0(h 2 ) if the coefficients 
of those terms, which involve V 4 <JG)>, are uniformly bounded. The biharmonic 
operator appears in the lowest truncation terms from Sn and the trapezoidal 

r' -sum. If Simpson’s rule for the r' -integration were used, there would be 
sixth-order derivatives in the lowest terms. The sharpness of the initial con- 
dition makes it difficult to maintain small truncation errors unless h « <7 0 and 
be chosen sufficiently small to strongly overdamp these errors. 

(2) Besides having an essential singularity at zero time, <(G)> has sharp and 
nearly discontinuous behavior that persists for slightly longer times, so that 


*We actually require G™ >N =0. These terms appear in H as formulated for the infinite domain and 
hopefully make a negligible contribution for large enought r^ . 
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nearly -unbounded derivatives are not restricted to t « 0, The Fourier trans- 
form of <G> for times t « \ / |u| is identical to the large Pecle't number result 
of Lee< 15, ff * EqUfltion < 29 >> for wavenumbers in the equilibrium subrange. In- 
stead of falling off monotonioally with increasing r , <(g)> rises to a cusp at r *= 

2 |u|t which becomes an infinitely sharp -cut-off when k *- Ch* 


lim 


<G (r, t)> s 


[277 |u| t]” 2 (4 |u | 2 t 2 - r 2 )" 1/2 , 
0 


r < 2 |u| t 
r > 2 |u| t 



« k ’ |u 


(3) The severity of the discretization error is sometimes illustrated by the 
appearance of large oscillations associated with the truncation of high order 
derivatives. Lack of oscillations does not imply convergence, but the appear- 
ance of large ones indicates that the discretization error is also large. In order 
to gain some insight into this behavior, consider the Crank-Nicolson case with 
uniform time increments. We assume that the nonhomogene ous part of Equation 
(1) can be expanded as « 2 s^X* and that G itself can be expanded as 
G n = ovor eigenvalues k j of the equation (A - I\j)X j = 0 where A is 

the matrix derived from the heat operator (with X N ~ 0) 


-6 

X 

2 ~2 

'X 0 

-12 2 

I 

CO 

to 

/ 

X *x 

s X 

X \ 

•XxX;0X 

X -( 

X 

*x. X. 

0 

x s x 'x s 
X ^(‘ -nTt) 


and the X’s correspond to the space-part of the separable solution to the heat 
equation. The problem is to determine the properties of affecting stability 


*This behavior of was discussed by Roberts^ 1 ^, who obtained the zero-K formula, and by 
Kraichnan^. The form of for k 0 corresponding to Roberts’ result has not been reduced 
beyond the integral forms ^G(r,t)^ =(277 2 r [u (t)“ Jq 03 dk sin(kr)J^ (2 |u |kt) exp{~rk 2 t) correspond- 
ing to direct inversion of Lee’s result, or 

<G(r,t)> s 

I (47rr)(2 | u 1 77 1 ) 2 V77*t 3" 1 u I 1 ' ) dt, [(2 |u| t) 2 .-£] h expfX^KtJs inh \x/Z 2* 1 1 

from a convolution theorem. 
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and oscillatory behavior. Straightforward substitution of these definitions into 
Equation (2) yields 




/ a j)m-m' yj 
nr 


where 


vj r . 
*0 


al ' K>o= + "V> 


a J 


1 - A-j p/( 1 + 0 k. p ) , p ~ K't/ h 2 , 


and a* are the coefficients defined by the initial condition 

N- l 

a* = r ^ G° X 1 

/ n n n 
n * 0 


f 


where are the weighting factors 


N- 1 



llaO 




K K 


a jk 


Crandall 1 * 1 7 * investigated the homogeneous system (H 0) for which T l = 

(a j ) m . For this special case stability obtains if | a | <1, and no oscillatory 
modes appear if a > 0; for given p , stable and non-osciilatory response ob- 
tains if 0 > 1 - l/p || K whereas mere stability requires 0 > 1/2 - l/p\\ \ .* 

The nonhomogeneous system treated here is far more complicated. The S*s are 
functionals of all the T’s and X,’s , and soT obeys a nonlinear functional equation 
with perhaps no uniform criteria for stability or suppression of oscillations. 

This analysis should be extended further and be applied to the multilayer and 
variable-^ forms of the difference equations. 

(4) The last expected difficulty that we mention concerns the rate of con- 

» 

vergence of the iteration sequence for the Crank-Nicolson form. Rapid con- 
vergence requires that z be a small fraction of the term €h ; but since H is a 
time convolution, sc is a linear functional of G° whicn is nearly singular. Con- 
sequently, z may make a significant contribution to and large initial 
truncation errors may propagate directly to g . 


* The uniform norm ||x|| TO •' maxj | \.| is called the spectral radius of A. Application of Gerschgorin’ s 
theorem yields an estimate of 17/2 as an upper bound of ||A|| (r) . 
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d, Summary of numerical work. 
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An incomplete and inconclusive investigation of numerical procedures for 
calculating G was carried out on the IBM 7040-7094 Direct-Couple System at 
the University of Washington Computer Center at various times during 1900 and 
19(57. Nearly all numerical experiments were made with the parametric values 
k a=r 1 and h » 0.1. The overall mesh size was (N, M, J) = (100, 200, 100) 
or (50, 100, 50) with nearly identical results in either case. 

Before the numerical schemes described above were tried, a straight- 
forward explicit scheme (0 =0) was tried with uniform •£ = 0.004 so that p = 2/5, 
the lower stability limit of the difference scheme used here for the heat equation. 
For | u 0 a strong dependence of the shape of G„ on the behavior of G was 
found. For | u | = 1 the scheme was extremely unstable, so this approach was 
abandoned. 

Next WPS tried the implicit scheme corresponding to Equation (2) with 
6 - 1/2, <w= 1.3, uniform & , and Gaussian G° with cr, = h . Fairly extensive tests 
were made to determine the accuracy of the numerical methods : one result was 
that the triple integration was only accurate to about 5%, regardless of whether 
the trapezoidal or Simpson’s rules were used. For | u | =0 and = 0.1 oscil- 
lations occurred with amplitudes and decay rates of the same order as for G 
itself; as £ was decreased, the amplitudes decreased until, for ■i - 0.01, G re- 
mained positive, and for 't = 0.004, the oscillations effectively disappeared. The 
calculation time for | u | / 0 was found to be excessive. For example, for the 
larger mesh it required 3.6 minutes to step off one time layer; whereas, for a 
course mesh with N = J = 25, h = 0.2, % = 0.02, and | u | =10, only 16 layers were 
calculated in 45 minutes using f 2 (r,t). No further experiments were made after 
it was found that H had been incorrectly reduced to the isotropic case for the 
implicit numerical method. 

It was decided to extend the computer program to include the variable 
and multilayer schemes before correcting H, since it was felt necessary to have 
small time increments initially and desxrable to eliminate the iteration. The 
main features of these programs are described in the Appendix. 


4 . SUGGESTIONS FOR COMPLETING THE CALCULATION 

The following suggestions may be helpful in reducing the severity of the 
numerical ti-oubles. 

(1) The dimensionality of the initial singularity may be too large. Lower- 
order delta-function initial conditions may be achieved, for example, by treating 
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r <b(r ,t)]> as the dependent variable. In this case the spherical Laplacian op- 
erator goes over to the one -dimensional Laplacian. 

(2) Because the initial singularity occurs at the origin r = 0, it may be 
appropriate to use nonuniform space increments hjj or expanded coordinates 
which would be time dependent to account for the spread of <£)>. The numerical 
analysis would be considerably more complicated than for uniform h . 

(3) It may be worthwhile to use the short-time analytic results for <jC)> to 
begin the calculation without spurious oscillations. The finite-difference calcu- 
lation would begin at some time t « k/ |u| where <^G)> is likely to be much 
smoother than it was initially. 

(4) The calculation can be done much more easily in Fourier coordinates 
for which the initial condition of <^G(k, t'p * is unity for all wavenumbers. The 
number of independent variables is the same for <Jg)> and<p^>, and both H and 
H require three integrations. In fact, Lee< 15 > has evaluated <(3 for f (r,t) 
corresponding to the energy spectrum function for the final period of decay and 
the asymptotic modal time correlation of Kraichnan/ 1 * We want, however, to be 
able to generalize the techniques to nonhomogeneous turbulence. For that case 
<(S)>is initially diagonal with respect to the two wavenumbers for the same 
directions of nonhomogeneity. Just as for the homogeneous case, and<^G)> 
have the same number of independent variables; H and H each have one time 
integration and the same number of integrations over coordinates or wave- 
numbers in directions of homogeneity. Unfortunately, for each nonhomogeneous 
direction, for which there is one integration in H, there are seven integrations (or 
sums for a finite domain) for H.^ For a suitable Fourier expansion of <uu') 

this seven-fold integration can be reduced to a three-fold one. If this number of 
integrations remains intolerable, an alternative procedure would be to calculate 
<(2^>for small times, invert to<(G^>, and then finish the calculation of <(G^> for 
intermediate and large times. 


*The tilde denotes a suitably-defined Fourier coefficient. 

^This can be seen by examining H, written symbolically as H = [V ‘(uu ><G> * [-V<G>] 
where the asterisk represents space-time convolution. For each nonhomogeneous direction H 
produces one integration from the space convolution, one from each gradient operation, and four 
from the <(uu ')> product. 
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APPENDIX 

a. Tabulation of truncation errors 

Direct application of the difference and quadrature operators in Section (3a) 
yields the following truncation errors: 


s n MO - MO = v ? n MO + 


1, n > 0 
^5 n 
3’ " = 0 


S n MO -37- MO = ^7 ^ MO + ••• . y “ 


X , n > 0 
11 = 0 


£ (O - MO h 2 JL-MO + • • ■ , i = x- 


mod x 


® r x 






<3) MO - MO =-| (O 2 — ^-MO + " ' . 

3 3t2 


A( t 2) «‘. fS >-3T-^( *-♦*> = *<*-'> 

m+9 3t m+0 

<U3£+3f) * )2 _if_ 0(t ) + ... 

(j v «n + l/ 3 m + 0/ T 

m + # 


A t 3) MO - af-MO 

dt m 0 B t 3 


and 


nm 


2^ P < r n’ r -’ O O V l<G»-H(r n , tj= © + © + © + © 


n ' j m* 
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where 


* 


r r N a 

® =h 2 2rr r * 2 dr' dfi dt' <^G(r', t')> 

4) J-l J 0 

+ T^t}< g( V *■-*'» ’ r y £ <°’ r N> 


© = 


7rh S 

18 




Jl 


P ( r n ’ f 1 ' r N^ 


® =(2/J) 2 


j>-j 

•/<> •'0 


P(r n , r‘, 77, t m , t' | <C>) , 7J£(-1, 1) 

OT) 2 


( 3m 2 \ z.® r l 

^ r ' 2dr ' L d/i ^ P(r " ,r '^’ t "- 5l < G »- wm.) 


Note that term (T) isO(h 2 ), © is 0(h 2 ) for finite r N> (§) is 0( [A/i] 2 ) which may 
be madeO(h 2 ) by taking J ~ 2N. The order of (T) depends onb and m. Term 
© isO(£ 2 ) as b - 1 andO([l + 3(b-l)' 2 } -1 ) as m - co for finite t m . The an- 
ticipation that 3 2 P/3t' 2 becomes small for t # near t m2 is the justification for 
assuming that the variable order of © can be effectively maintained atO(£ 2 ). 
With these restrictions the truncation error for H is 0(h 2 ) + 0(£ 2 ). 

b. Thomas algorithm 

Consider Equation (4) with z lumped into d; i.e., Bg = d. B is tri-diagonal 
with diagonal elements b 0 , , . . . b N-1 , subdiagonals -a 1 , - a 2 , . . . - a N _ t , and 

supradiagonals - c 0 , - c t , . . . - c N _ 2 . Normalized upper-diagonal form is obtained 
by calculating new coefficients 


°0 = c 0 /b 0- % = c n /e n 


4> = V b 0' = ( d „ + a n ^n-l) /e „ 


r n = 1 , 2 , 


J 


N - 1 
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where e n = b n - a n a n _ r Then g is given directly by 


g N = °> gn = B„ + l °n + / 3 n n = N - 1, N - 2 , • • • 1.0. 


c. Computer programs 

The computer language used for the numerical experiments is Fortran -IV, 
Because of storage limitations (core = 2 15 = 32,768 words) the calculation of G 
is performed in one job segment in ALTIO mode (minimal 10 buffer), and the 
result is written on magnetic tape; in the second job segment G is read from 
the tape in FIOCS mode and is printed, plotted, or used in subsequent calcula- 
tions, To avoid long subroutine argument lists, G , all input parameters, and 
certain calculated quantities are placed in a COMMON block; those portions 
which are used in each subroutine are so indicated by EQUIVALENCE state- 
ments. In the first job segment G is stored in a 102 x 201 array. 

In the first job segment the magnetic tape is first positioned with a tape 
manipulation subroutine. Next a subroutine for handling input reads in the 
necessary input parameters mentioned in the text, computes auxiliary param- 
eters, and places them all in COMMON. This routine also has an entry pro- 
vision for changing input parameters under a NAMELIST format. Then a sub- 
routine is called which performs the computation of G using either Equations (2) 
or (3) and the methods of Section (3b). This routine calls routines for calculating 
Gn and the term (d + z ) in Equation (4). 

In the second job segment the values of G calculated in the first segment are 
either printed out or jotted (with the University of Michigan plotting routine) or 
else made available for further computations. 
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